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SUMMARY 


A systematic procedure is developed for the calculation of the 
structural response of aircraft flying through a gust by use of differ- 
ence equations and matrix notation. The use of difference equations in 
the solution of dynamic problems is first illustrated by means of a 
simple -damped-oscillator example. A detailed analysis is then given 
which leads to a recurrence matrix equation for the determination of 
the response of an airplane in a gust. The method takes into account 
wing bending and twisting deformations, fuselage deflection, vertical 
and pitching motion of the airplane, and some tail forces. The method 
is based on aerodynamic strip theory, but compressibility and three- 
dimensional aerodynamic effects can be taken into account approximately 
by means of over-all corrections. Either a sharp-edge gust or a gust 
of arbitrary shape in the spanwise or flight directions may be treated. 
In order to aid in the application of the method to any specific, case, 
a suggested computational procedure is included. 

The possibilities of applying the method to a variety of transient 
aircraft problems > such as landing, are brought out. A brief review of 
matrix algebra, covering the extent to which it- is used in the analysis, 
is also included. 

INTRODUCTION 


In the problem of an airplane flying through gusts, accurate 
predictions of stresses are not always obtained if the interaction 
between aerodynamic loads and structural deformations is not considered. 
The present paper gives a method for determining the dynamic response 
of aircraft in gusts in which this interaction is considered. An 
approach is employed which is a departure from the usual modal type of 
solution. The time derivatives in the integro-dif ferential equations 
of motion of the airplane are replaced by appropriate difference 
expressions and use is made of matrix notation to express conveniently 
the conditions of equilibrium at a number of points along the wing span. 
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The result is a systematic procedure which is complete and general in 
form. The airplane is assumed to be free to translate and pitch. Wing 
bending, wing twist, and fuselage flexibility are all included. Tail 
forces due to vertical motion, angle of attack, and gust penetration are 
also included in the analysis. 

With the method, a gust with any gradient in the direction of 
flight or along the span may be treated without difficulty. The method 
is based on aerodynamic strip theory, but over-all compressibility and 
aspect -ratio corrections may be included without difficulty, if desired. 
One such over -all correction is indicated. 

In the first part of the paper the method of using difference 
equations in the solution of dynamic problems is illustrated by an 
example in which the transient response of a simple oscillator is 
determined. The analysis for the determination of the response of an . 
airplane in a gust is then given. In the following section a computa- 
tional procedure is suggested. This section is not intended to 
describe or add to the understanding of the analysis, but by following 
the directions indicated, the response of any airplane may be found 
without following through the complete details of the analysis. 

Supplementary definitions and derivations are presented in appen- 
dixes. Appendix A summarizes the various matrix coefficients and 
matrices that are used in the analysis, appendix B gives a derivation 
of the difference equations, appendix C gives a derivation of the 
flexibility matrices, appendix D gives a derivation of a recurrence 
equation for evaluating the form of Duhamel's integral which involves 
an exponential kernel, and appendix E presents a review of the funda- 
mentals of matrix algebra. It is suggested that those not familiar 
with matrix algebra read appendix E before reading the analysis. 


SYMBOLS 


a distance between leading edge of wing and elastic axis 

8q_ coefficient used in unsteady lift function for sudden change 

in angle of attack 

A aspect ratio of wing 

Aj. aspect ratio of horizontal tail 

b semi span of wing 

c chord of wing 
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C °t 

c t 


E 

F 

G 

i 


I 

J 

k 

Z 

L 



chord at wing midspan 
midspan chord of tail 
mean aerodynamic chord of tail 

distance between mass center of wing cross section and 
elastic axis of wing; positive when elastic axis lies 
forward of mass center 

Young's modulus of elasticity 

suddenly applied force 

shear modulus of elasticity 

integers 0, 1, 2, 3, and 5 used to designate stations 

(for most part used as parenthetical numbers, that is, w(3) 
is deflection at station 3) 

bending moment of inertia 

torsional stiffness constant 

radius of gyration of wing mass about elastic axis or 
elastic spring constant 

length of section associated with a spanwise station 

aerodynamic lift over interval Z on wing 

shear force transmitted to wing by fuselage 

aerodynamic lift over interval Z on wing due to gust 

one -half aerodynamic lift on tail due to gust 

one -half total aerodynamic lift on tail 

part of aerodynamic lift over interval Z on wing (see 
equation (l6)) 

part of aerodynamic lift over interval Z on wing (see 
equation ( 17 )) 

mass of beam included in interval Z or concentrated mass 
in spring oscillator 


m 



k 
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m 


m A 


m At 


me 




M 

M f 

n 

P 

P f 

q 

s 

As 

St 

t, T 


/ 2 \ 

mass' m including apparent mass effect f m + j 

assumed over-all compressibility and aspect-ratio correction 
A 


for wing 


(2 + A\jl - M 2 ) 


assumed over-all compressibility and aspect-ratio correction 

for horizontal tail ( f \ 

\2 + A+Wl - M 2 / 


mass moment me including apparent mass 


effect me + 


(- + - 1 )) 


mass of fuselage per unit, length 

mass polar moment of inertia mk 2 including apparent mass 

' M j. / ,2 JtpZcVl a\2 npZc^ 

effects (mk + - ~) * -^§5 

Mach number or aerodynamic moment over interval l about 
elastic axis of wing 

moment transmitted to wing by fuselage 

integers 0 , 1 , 2, 3 , and so forth to designate number of time 
intervals passed 


normal load acting at a station 
fuselage inertia loading per unit length 
torsional load acting at a station 
distance traveled by wing in half -chords 


chord c Q is used as reference chord j 


(?*■ 


where midspan 


(SO 


distance interval in half -chords 
horizontal -tail area 

time, zero at beginning of gust penetration 


/ 
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U forward velocity of flight 

v vertical velocity of gust 

w deflection of elastic axis of wing, positive upward, or 

deflection of mass oscillator 


Wf 

W X 


x 


'•n 


x t 


deflection of fuselage, positive upward 

fuselage modal function, zero at wing elastic axis and unity 
at tail one -quarter -chord location 

distance along fuselage measured from wing elastic axis , 
positive in rearward direction 

distance from foremost part of nose to elastic axis 

distance from elastic axis to one -quarter -chord location on 
tail 


y 

z 

a t 

0 

0t 

7 

6 


distance along wing measured from center of airplane 
ratio of dynamic deflection to static deflection 
angle of attack of horizontal tail 

forward-speed and aspect-ratio factor for wing (m^rtpU) or 
coefficient of damping for spring oscillator 

forward-speed and aspect-ratio factor for tail (h m^npS-tU 
exponential coefficient in 0 function associated with 

tuna t ,(7 = f o x) 

coefficient of fuselage modal function 


€ time interval 

\ exponential coefficient in 0 function associated with 

variable s 

dimensionless interval between i - 1 and i 
stations (Xjb is actual length) 

p mass density of air 

<p . angle of twist of wing, positive in stalling direction 
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\|r function which denotes growth of lift on rigid airfoil 

entering sharp-edge gust (used without subscript to 
indicate function for wing and with subscript t used 

to indicate function for tail) 

/ 

u>f natural frequency associated with Wq, radians per second 


lit) 

1-4 


[] 
rj 
1 ! 

LJ 


unit -step function 

function which denotes growth of lift on airfoil following 
sudden change in angle of attack (used without subscript 
to indicate function for wing and with subscript t used 
to indicate function for tail) 

square matrix 

rectangular matrix 

column matrix 

row matrix 


Subscripts: 

t tail 

0,1, 2, 3, . . . n number of time intervals passed 

0,l,2,3,k,5 or i station (however, station is usually given as 

parenthetical number, such as w(3) for deflection 
at station 3)} i is also used as general subscript 
in appendix A 


All the terms, coefficients, and matrices not defined in this 
section are defined in appendix A. 

Dots are used to indicate derivatives with respect to timej for 

5w • Sw * 

example, ^ = w or ^ =? w . 

ANALYSIS 

Transient Response of a Simple Damped Oscillator 


In order to illustrate the use of difference equations and to test 
the accuracy, of the procedure that is to be used in the solution of the 
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more complicated gust problems, the solution of a simple problem having 
a known analytical solution is first presented. The problem is to com- 
pute the response of the damped oscillator shown in figure 1 to a 
suddenly applied force. The differential equation of motion of this 
system due to the suddenly applied force is 

mw + (3w + kw = F/( t) (l) 


By use of difference equations this differential equation may be trans- 
formed into an equation which involves deflection ordinates at several 
successive values of time. Probably the most commonly used difference 
equations are the following (see appendix B for derivation): 



v n+l ~ v n-l 
2e 


( 2 ) 


w n = 


w n+l ’ 2w n + w n-l 


(3) 


which give the derivatives at the intermediate of three successive 
ordinates. Although these equations are quite adequate for the 
oscillator problem of the present paper, they cannot be used in the 
gust analysis which follows. Rather, for reasons which are brought out 
in a subsequent part of the analysis, equations that give the derivatives 
at the end ordinate of several successive ordinates must be used. If 
only three successive ordinates are used, the derivatives so found are 
not accurate enough to be useful. If a fourth ordinate is added, how- 
ever, derivatives may be taken at the end ordinate with accuracies 
which are comparable to those given by equations (2) and (3). Such 
derivatives are derived also in appendix B and are given by the 
equations: . . 


"n 


llw n - 18*^ + 9w n _ 2 - 2w 3 

■ 


w 


Wn = 


2w n - gWn-i + W n _ 2 - w n _- 


(5) 


Although either equations (2) and (3) or equations (4) and (5) may be 
used in the solution of this oscillator problem, only equations (4) 
and (5) will be used, since only these equations can be used in the 
gust-problem solution presented in this paper. 
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If the derivatives in equation (l) are replaced by the difference 
equations (4) and (5), the following equation is obtained: 



6 m 





m 


v n-l 





( 6 ) 


This equation may be said to be the difference equation of motion. If 
the more general case of a variable applied force were- being considered, 
the factor F in this equation would be replaced by F n , the value of 
the force present at the time t = ng. 

k B 

If a specific case is now considered, in which — = 400, — = 2, 

€ = 0.01, F = 1, and the notation z =' (ratio of dynamic deflection 
to static deflection) is used, equation (6) becomes 


z n = 0.018927 + 2.42272 z n -i - 1.92114 z n _2 + 0.47949 z n _3 (7) 

This equation may be regarded as a recurrence formula} the value z n 
may be interpreted as the deflection to come and may be found easily 
from the three preceding deflections z n _i, Zn_2.> and Zn_ 3 . Then with 
the newly found value z n and with z n _^ and z n _ 2 , the next deflec- 
tion can be found, and so on. This process thus gives a step-by-step 
derivation of the time history of deflection and may be carried out as 
far as is desired. Of course the process must start with known initial 
values of z. These values can be found with the aid of the initial 
conditions of the problem by means of the following approach. 

The difference equations for the first and second derivatives at 
the third ordinate of four successive ordinat'es are (see appendix B) 


wh = -^(2w n +l + 3wn - 6wn-l + Yn-2^) 

"n = -^(v n+1 - 2w n + w n _q) 

If these equations are taken to apply at t = 0 (n = 0), they become 
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wq = £“(^1 + 3wq - 6w-i + w_ 2 ^) (8) 

% = -^( w l " 2w 0 + w_J (9) 

6 

For the problem under consideration the primary initial conditions are 
that, at t = 0, the displacement and velocity are zero. By. use of 
equation (l) or by reasoning from Newton's second law, a secondary 
initial condition can be established* that is, the acceleration 
immediately following the application of the unit force must be l/m. 

In equation form these conditions axe 

w 0 = 0 

w 0 = 0 



By substitution of these values into equations (8) and (9) and by use 
of the notation z = the following relations can be found to exist, 

between the ordinates: 

Z 0 = 0 


z -2 = 0«24 “ 8z x > (10) 

Z -1 = 0*0^ " z i 

J 


Substitution of these values into equation (7)> with n set equal to 1, 
gives an equation from which zq, the deflection at t = € , may be 
evaluated. Three successive deflections can now be established: the 

deflection at t = e , the zero deflection at t = 0, and a fictitious 
deflection for t = -€ given by equation (10). In the real problem no 
deflection exists for t less than zero; the assumption that^ a deflec- 
tion does exist before t is zero is simply a means for allowing the 
recurrence formula, equation (7)> to apply at the origin as well as at 
later values of time. Furthermore, no violation is made of the condi- 
tions under consideration because, mathematically, the response 
after t = 0 is not influenced by the response that may exist 
before t = 0, so long as the initial conditions are satisfied. The 
process just described for treating the initial conditions is actually 
not different from the . procedure commonly employed in difference- 
equation approaches, in which exterior points near a region under con- 
sideration are written in terms of the interior points by means of the 
boundary conditions. 
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With the initial values of deflection thus established the step- 
by-step evaluation of succeeding deflections proceeds in a straight- 
forward manner; that is, equation (7) is now evaluated for n = 2, 
then for n = 3, and so on. The response of the oscillator for the 
physical constants listed previously is given in figure 2. The compari- 
son between the difference solution shown in this figure and the exact 
solution of equation (l) is seen to be good. As a matter of interest, 
the solution is also shown in this figure that is obtained by the use 
of the parabolic end-ordinate derivative which involves only three 
successive ordinates. The agreement in this case is seen to be quite 
bad. If equations (2) and ( 3 ). had been used, on' the other hand, the 
difference solution (in this case for w n+ ^) would correspond to that 

given for the cubic end -ordinate derivative. 


Recurrence Matrix Equation for Response of 
an Airplane in a Gust 

In order to help the reader to obtain a perspective of what is to 
be covered in this section, the following basic phases .of the analysis 
are given: 

(1) The assumptions are stated. 

(2) The coordinate system and displacements are defined. 

(3) The aerodynamic lift and moment are defined. 

(4) The normal and torsional dynamic loadings (inertia forces, 
aerodynamic forces, and fuselage forces) on the wing are derived. 

(5) The equations of elastic deformation - wing vertical motion, 
wing rotation, and fuselage bending - are given. 

(6) The dynamic loadings on the wing are transformed into 
difference equations. 

(7) The equations of elastic deformation and the difference 
equations for loading are combined to give the recurrence matrix 
equation for response. 

In succeeding sections the initial response is determined, the 
method for evaluating the gust forces is shown, and the method for 
computing the loads and stresses is indicated. 
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Assumptions . - In this analysis an attempt is made to obtain the 
simplest and most direct solution to the problem with a minimum of 
simplifying assumptions. The case treated is that of an airplane having 
an essentially straight wing and penetrating a gust of known structure. 
The tail is considered to penetrate subsequently the same gust as does 
the wing. The assumptions made are as follows: 

Assumptions pertaining to elasticity and airplane motion: 

(1) The usual assumptions of engineering beam theory are made. 

(2) The fuselage is free to pitch and move vertically. The portion 
of the fuselage in front of the elastic axis of the wing is assumed for 
convenience to be rigid. The portion of the fuselage rearward of the 
elastic axis is assumed flexible, and the elastic deflection is assumed 
to be given by a single modal function. 

(3) The tail is assumed rigid. 

Assumption's pertaining to aerodynamic forces: 

(1) Aerodynamic strip theory applies. Three-dimensional effects, 
however, may be taken into account approximately by means of over-all 
corrections. Some such corrections are indicated. 

(2) The gust force and forces due to vertical and pitching motion 
are the only tail forces considered. Other forces of known c har acter 
may be included, however, if desired. 

(3) Aerodynamic, forces on the fuselage are neglected. 

Coordinate system and displacements .- Position on the airplane is 
denoted by an orthogonal system of axes. The origin is at the inter- 
section of the wing elastic axis with the plane of symmetry of the 
airplane: the w-axis runs positive upward, the x-axis runs along the 

fuselage positive in the rearward direction, and the y-axis runs span- 
wise. The wing semispan is considered' to be divided into six, not 
necessarily equal, sections, with a station point at the middle of each 
section. (See fig. 3.) More or fewer stations could be chosen, but it 
is believed that six is a fair compromise between the amount of labor 
involved in setting up a solution and the accuracy desired. The interval 
between stations is designated by the number of the station at the out- 
board end of the interval. Station 0 is located near the wing root and 
generally may be located where the fuselage intersects the wing. In this 
way the concentrated forces of the fuselage are allowed to act through 
station 0. The other five stations are then located in any convenient 
manner so as to fall at concentrated mass locations or at points which 
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represent the average of distributed masses, station 5 being nearest 
the tip. The total mass within a section is assumed to be concentrated 
at the station point, and the average of the section geometry (chord, 
elastic axis position, and so on) is assumed to apply. In this way the 
wing is assumed to be a beam subject to six load concentrations, and as 
such will have a linear moment variation between each station. The 

further assumption is made that the variation is linear between 

El 

each station. With these assumptions for the El variation and con- 
centrated load locations, equations for deflection at each station 
point may be derived (appendix C) by direct analytical treatment. 

The displacements of the cross section at each station of the wing 
are given as the deflection of and rotation about the wing elastic axis 
as shown in figure 4. The fuselage displacements are shown in figure 5 
and are given by the equations: 

p 

' w f = w(0) - cp(0)x (11) 

for the forward section and 

Wf = w(0) - cp(0)x. + W}_& (12) 


for the rearward section. The function W]_ is taken as the fundamental 
mode of vibration of the fuselage and tail assembly, when the fuselage 
is considered to be clamped as a cantilever beam at the elastic-axis 
location of the wing, and is given in terms of a unit deflection at 

the ^-chord position on the tail. With this function to represent the 

elastic deformation of the fuselage the deflection and angle of attack - 
of the tail is found with the aid of equation ( 12) to be 

w f (x t ) = w(0) - qp(0)x t + 6 ( 13 ) 


where 


«t = - 


dWf-1 


A 



= <p(0) - 80! 

J 


dW x 
Q 1 = dx - 


( 14 ) 


x=x t 
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Aerodynamic lift and moment . - Before . going into the details of the 
analysis it is felt worthwhile to define and describe the nature of the 
aerodynamic forces to which the wing is subjected. These forces 
originate from two sources: they arise directly from the gust encoun- 

tered, and they arise from the ensuing airplane motion. The equations 
for the aerodynamic lift and moment that develops are herein set up in 
a convenient form on the basis of work given in references 1 to 4. In 
these investigations various methods for separating the lift forces 
have been used, but the particular method for separating these forces 
is not important so long as they are taken into account properly. 


In the present paper the aerodynamic lift and moment are considered 
to be composed of two parts: one part, designated as the lift or 

moment due to circulation, which includes aid. lift forces or moments 
exclusive of aerodynamic inertia effects and the other part, which is 
due solely to these inertia effects. These lift forces and moments can 
be resolved into the force systems acting on the airfoil as shown in 
the following sketches: 


Forces due to circulation 


Inertia force and moment 



The force Lg is the lift force developed by the gust. All the other 
forces occur as a result of motion of the airfoil. These forces, as 
well as the gust force, are given for an interval Z of the span by 
the equations: For the forces due to circulation, . 


r»t- ^ 

L g = mA«pcZU J \|r(t - t ) dT 


(15) 


r* r. •• fa a V- 

1 = m A npcZU J U<p - w + c^f - ~ycp 


1 - <J>(t - t) 


dT 


(16) 
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lk 


L2 = - 


m A rtp 2 c c 


Ucp 


(17) 


and for the inertia force and moment. 




-w + I — 


a) 9 


(18) 


Ml = ‘ 2 ®" * (19) 


factor which can be used to introduce over-all compressibility 
and aspect-ratio corrections} in this paper the factor is 

A 

assumed, to be given by ■■ ■ 

2 + A\jl - M 2 

1 - $ lift function which denotes the growth of lift on an airfoil 
following a sudden change in angle of attack 

t lift function which denotes the. growth of lift on a rigid 

airfoil entering a sharp-edge gust 

The functions 1 - S> and ijr and the correction m A are estab- 
lished as follows. In reference 5, approximate equations are derived 
which give the lift -coefficient form of the growth of lift on a finite 
wing following a sudden change in angle of attack or due to the penetra- 
tion of a sharp-edge gust. The equations may conveniently be considered 
as the product of a factor, which may be regarded as a lift-curve slope, 
and an unsteady lift function, designated by 1 - $ for the function 
due to the angle rof -attack change and by \|/ for the function due to 
the sharp-edge gust. These unsteady lift functions are shown in 
figures 6 and 7 and are given by the following equations: For 

the 1 - $ functions 


(1 - $) A=3 = 1 - 0.283e" 0,540s (20a) 
(l - <5) a= £ = 1 - 0.36le"°* 38ls (20b) 
(1 - 4>) a=00 - 1 - 0.l65e'°*° 45s -0.335e"°’ 3 ° OS 


where 


m. 


(20c) 
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and for the \|r functions 

'♦a= 3 = 1 - 0.679e"°* 558s - 

* A= g = 1 - 0 .448e"° * 290s - 
\ =OB = 1 - 0.236e ^ - 

\ = 1 - 0.500e"° ,13 ° S - 

A= 00 


-3.20s 

0.227e 


(21a) 

0 .272e’° *^ 25S - 

0.193e‘ 3 - 00s 

(21b) 

0.513e'°* 364e .- 

0.171e' 2 ' , * 2s 

(21c) 

0*500e " S 


(22) 


Equations (21) are based on equations of reference 5) whereas equa- 
tion (22) is the i|r function that is suggested for wings of infinite 
aspect ratio in reference 3- Inspection of equations (20) shows that 
the $ function for aspect ratios 3 and 6 is given by a single 
exponential term. It is probable that the 0 function for higher 
aspect ratios, say 10 and even 20, may also be given to a sufficient 
approximation by a single exponential term. Therefore, the assumption 
is made that in general .G> may be represented by an equation of the 
form 


® = a ie " Xs . (23) 

Interpolation, for example, of the curves in figure 6 shows that 
the $ function for aspect ratio 10 might be approximated by the 
equation: 

\=10 = °- 4le "°' 3S (24) 

The analysis does not necessarily limit $ to a single exponential 
term. Other terms could be added with some increase in labor, but the 
degree of refinement obtained is not expected to add much to the over- 
all accuracy of the solution. 

Although the functions given by equations (20) to (22) are known 
to approximate the true functions quite well over a large range 

tj , the \|/ functions given by equation (21) do not 

vanish, as they should, when t = 0. When used in the computational 
procedures given hereinafter, these functions, therefore, are to be 
taken as zero when t = 0. Another point to note is that the variable s 
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is given in terms of a reference chord c Q ; thus this variable as 
applied to the wing is different, in general, from the variable as applied 
to the tail. 

Examination of the values of lift-curve slope, which were stated 
to be present in the equations taken from reference reveals that they 
may be approximated with good accuracy by the product of 2 it and the 

A 

often-used aspect-ratio correction — — — for steady incompressible 

flow. In the present paper it is assumed that compressibility and * 
aspect-ratio corrections can be made by replacing this aspect-ratio 

A 1 

correction by a compressible aspect-ratio correction defined by — — 

where A' = A\jl - M^, and by multiplying this -correction by the Glauert- 

Prandtl Mach number correction . ■ to give the product m^. The 

V 1 - M 2 

procedure then for taking into account three-dimensional and compressi- 
bility effects in’ the present analysis is to determine m^ from the 
forward speed and aspect ratio of the wing and to use the 1-0 
and r functions, equations (20) to (24), for the aspect ratio which 
is nearest that of the wing. 


Some word of explanation of equation (l6) might be worthwhile at 
this point. The 0 (t - t) function is associated with the lift forces 
which are due to the wake. Without this term the equation would yield 
the steady lift corresponding to the instantaneous values of angle of 
attack and vertical velocity. If equation (l6) is integrated by parts 
and the conditions are stipulated that w, w, cp, and cp are all zero 
at t = 0, the following equation may be found: 


|3cZ4qW - (l - 4q)0cZw + PcZU 1-4 


E 



(1 - $q)0c^Z^ - ^ip + 3cZ J' vi(t - r)dr - $clU cpi(t - T)dT 


ec 2 z(3 - 



- t)1t 


(25) 


where (3 has been introduced as a forward-speed and aspect-ratio 
parameter defined by the equation 


(3 = m^TtpU 


(26) 
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With reference toequation (23 ) > $0 ^0 ■*’ n equation (25) would 

have the values 


$0 = 8.1 


2U , 

= * — ^1 
c o 


The form of Lq given by equation (25) is presented because this form 
is more, convenient to use in the present analysis. 


For this analysis the total lift and moment acting at the elastic - 
axis location are desired. For the present, the total lift L and 
moment M of the forces due to circulation are found; the inertia 
force and moment are to he treated separately. Summation of all the 
lift forces due to circulation and summation of the moments of *these 
forces about the elastic axis gives the desired equations for the aero- 
dynamic lift and moment acting on the airfoil over an interval l as 
follows: 


L = L x + L 2 + L g • (27) 

M = (a. - - (if " + ( a ‘ f) L g ( 28 ) 

The loading on the wing . - The normal and torsional dynamic loads 
that are held in equilibrium by the elastic restoring forces of the 
wing may be found by considering all the aerodynamic and inertia forces 
that act on the wing. The mass situated at any station (see fig. it-) 
can be shown to have an inertia normal force equal to 


-mw + mecp 

and an inertia torsional moment about the elastic axis equal to 

mew - mk^qj 

If the aerodynamic forces and moments (see equations (18), (19) , (27), 
and (28)) are added to these inertia loadings, the total normal and 
torsional loadings on the wing at each station are found to be given, 
respectively, by the equations: 

p = -mw + meq5 + L + 

q = mew - mk^qJ + M - - ajlL^ + Mq 
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The terms L3 and Mi ordinarily would appear with the aerodynamic 
lift and moment values but axe treated separately so that they can be 
combined with the structural mass terms. If use is made of equations (l8) 
and (19), the loading equations become 


p = - mv + me<p + L 

(29) 

• # “O . • 

q = mew - mk c cp + M 

. (30) 


where 



The terms appearing with the structural mass quantities in the defini- 
tions of .5, me, and mk 2 are the terms which are commonly associated 
with apparent mass effects. 

The value of the shear forces Lf and the moment Mf transmitted 
to the wing by the fuselage can be found in the following manner: From 
equations (11 ) and (12) the values of the inertia loading on the forward 
and rearward sections of the fuselage can be shown to be given, respec- 
tively, by the equations: 

p f = -m f p( 0 ) - cp( 0 )xj (31) 


Pf “ “IDf 



X + 



(32) 


Integration of these inertia loadings over the length of the fuselage 
and addition of the aerodynamic tail load 2Lt give the value of the 
total load' transmitted to the wingj one half of this load is designated 
by Lf and is assumed to act at station 0, the other half being con- 
sidered to act through the corresponding station on the other half of 
the wing. Integration of the moment of the inertia loading about the 
elastic-axis location and addition of the moment -2xfLt of the tail 
forces give the total moment due to the fuselage j one half of the 
moment is designated Mf and acts at station 0. The values of Lf 
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and Mf thus found can be given by the equations: 

L f = -M-lw(O) + MgCp(O) - M 38 + Lf. (33) 

Mf = M 2 w(0) - M4Cp(0) + M 56 - x t L t (34) 


where the Mf 's sire considered to be generalized masses defined as 
follows: 


1 l' x t 

a l = 2 


\ 


M i = — / mf dx 


1 r x t 

M 2 = 2 / m f x 

J *n 


1 r x t 

Mo = 2 / mfWf dx 
J 0 


M 4 = | r m f x 2 dx 

Jx n 


i r x t 

^ J m^xW^ dx 


1 r x t 2 

M 6 = | J Q mfWfdx 


V 


(35) 




The generalized mass constant Mg, although not appearing in equations (33) 
or (34), is included in this group because it occurs in a subsequent part of 
the analysis. In the derivation of equation (34), the aerodynamic moment 

of the tail about the tail ^-chord position is neglected since it is 

considered to be small in comparison with the value x^Lf.. The lift on 
the tail Lg can be found by application of equation (27) to the tail 
surface. In this case the <J> function appropriate to the tail should 
be chosen and the values of displacement w and q> should be replaced 
by Wf(xg) and at, the values of deflection and angle of attack at 
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the tail ^-chord position. These values are given by equations (13) 
and ( 14 ) . 

Matrix equation of equilibrium . - The problem of computing the 
response may be considered to be one of the determination of the deflec- 
tion and. rotation of a beam which is subjected to normal and torque t 
loadings. In differential form, the bending and rotational displace- 
ments are related to the normal and torque loadings by the well-known 
expressions: 


dy 2 


dy 2 


p 


(36) 



(37) 


where in this instance p and q are the loadings per unit length of 
beam. In addition to these two . equations which are considered to apply 
to the wing, an equation for computing the elastic deformations of the 
fuselage may be found} this equation may be found in the following manner 
The rearward part of the fuselage is considered to be a cantilever beam 
subjected to the inertia loading given by equation (32) and the tail- 
force 2Lt* If equation (36) is applied to the fuselage and use is made 
of equations (12) and (32), the following equation for fuselage bending 
results: 


5 


a 2 ' 

ax 2 EIf ax 2 



- cp(0)x + bW^ + 2L^. 


(38) 


in which Lt must be treated properly as a concentrated load and If 
is the bending moment of inertia of the fuselage. Since Wl represents 
a vibration modal function, the following relation exists: 


2 aV 

a^ 2 EIf a^ 2 ” 


= mfO)f Wf 


where o>f is the frequency of vibration associated with Wq. Equa- 
tion (38) may therefore be written 

brnfOi^W-L = -nif fi?(0) - q5(0)x + 5Wq”| + 2L t 
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Multiplication of this equation through by Wi and integration 
between 0 and xt re suit s in the following equation for fuselage 
bending 

oif^MgS = -M^w(O) + Mc>’(p(0) - Mg5 + (39) 

where M3, M5, and Mg are defined by equations (35) • 

Equations (36), (37), and (39), when the loadings given by equa- 
tions (29) and (30) are considered, are seen to be rather involved 
integro-differential equations but describe completely the motion of 
the airplane. The problem is to find functions w, 9, and 5 which 
satisfy these equations and which satisfy both the boundary conditions 
and the initial conditions. 


The problem of finding the w and cp functions may, be simplified 
considerably by reducing the rather complicated equations of motion to 
a simplified and systematic algebraic form. The first step (see 
appendix C) is to replace the differential equations (36)' and (37) for 
wing deflection and wing rotation by the following simple matrix 
equations: 


l_A]|w| = |p 

MM = h 


(bO) 

(bl) 


The matrices in these equations are defined in appendix C (see equa- 
tions (C22) and (C23) and equations (C29) and (C30), respectively) and 
have been derived on the basis that the displacements along the semispan 
are given at six stations. 


Equations (40) and C^l) and the fuselage deflection coefficient 5 
are now combined in a single matrix equation of the form indicated as 
follows : 



(te) 


This form is chosen because it will be useful subsequently. With the 
notation given in appendix A, equation (42) may be abbreviated to the 
form: 
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(43) 


This equation may be regarded as the loading matrix equation of equi- 
librium; it relates the loadings to the displacements by linear simul- 
taneous equations. The boundary conditions are automatically satisfied 
when this equation is used because they had to be taken into account 

when the submatrices [aJ and [bJ were derived. Only the initial 
conditions remain to be satisfied and these are treated separately in a 
subsequent section. 


Transformation of the loading equations into difference form . - 
The loading equations are now simplified by replacing the time deriva- 
tives by difference equations: If equation (5) is used to replace the 

derivative in equations (29) and (30) the values of the loading at 
the nth time interval are found to be 


Pn= “"7>(2w n - 5w n _! + 4w n _ 2 - w n _ 3 ) + ~(2^> n - 5q> n -l + 4qp n _ 2 - cp n _ 3 ) + I* 
€ € 

(44) 

^-(2<P n - 5q> n -l + 4<P n -2 " <Pn- 3 > + “n 

(45) 

The values L n and M n are found by determining the expressions 
for L]_, L 2 , and L g at t = ne (see equations (27) and (28)). Of 

these L]_ is the most complicated, since it (see equation ( 25 )) involves 

three unsteady lift integrals of the Duhamel type. Fortunately, however, 
,a rather simple recurrence relation can be developed which allows the 
calculation of the value of these integrals at a given t im e interval 
directly from the value at the previous time interval. This derivation 
is presented in appendix D and is made' possible because the <D function 
is of an exponential form. (Where the 0 function is given by more 
than one exponential term, a recurrence relation for each term may be 
written.) From the derivation in appendix D, therefore, the value of 
the three integrals at the nth time interval may be given as follows r 


q n= -^( 2v n “ 5w n _! + 4w n _ 2 - w n _ 3 ) - 


F 


n 


+ 


(3cZ e 
2 





(46) 
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where 


F = e 
r n c 


-7e 


F n-1 + e^n-l + §%-! 


in which g and g' are defined by equations (A5) in appendix A. With 
this expression to replace the value of the integrals in equation (25 ) } 


the value of Li n may be written 


ju(i - ® 0 ) - 4r b 0 - 


L^ n = PcZ^$o + § ®o) w n $o)P c ^n + 

(<j>0 + f $o) c (if ’ ?],n + f 1 " S’o) 1302 ^ " f)^n + F n ^) 

With the use of difference equation (4), this equation may be trans- 
formed finally into the form: 


L 1 = d 0 w n + d l w n-l + d 2 v n-2 + d 3 w n-3 + d 0% + d l%-l + 


n 


d 2 ,( Pn-2 + d 3 ,( Pn-3 + F n 


(48) 


where the d's are defined in appendix A. Likewise, from equations (4), 
(17), and (26), L2 n may be written 

L 2 n = §f-( 11( Pn - iaPn-1 + 99n-2 - aPn-3) 

If Li , L2 , and the value Lg n of the gust force at t = ne are 

introduced into equation (44), the value of p at the nth time interval 
can be shown to be given by the equation: 


p n " Vn * Vn-l* Vn-2 + Vn-3 + V’n + Wl + \ V 2 + 


l3%-3 + F ri + 1 


Sn 


(^9) 


where the tj's are coefficients which are given by equations (A3) in 
appendix A. In a similar manner, the equation for q (equation (^-5)) 
can be reduced to the form: 

* Vn + V l u n-1 + Vn-2 + Vn-3 + V^n * \ Vl + W 2 * 
v 3 '"Pn-3 * (■» - fK * (a - ^)tg n 

where the V*s are given by equations (A4) in appendix A. 


( 50 ) 
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The value of aerodynamic lift acting at the tail yj-- chord Lt is 

found most conveniently by applying equation (49) to one -half of the 
tail surface. This application is made by modifying the q coeffi- 
cients as follows: The mass value m is set equal to zero, — is 

1 c 
taken as c is replaced by c-fc, and 3cZ is replaced by defined- 

as the forward-speed and aspect-ratio parameter of the tail by the 
equation: 


h = \ m At «pUS t (51) 

In addition, w and cp are replaced by the deflection and rotation of 

the tail given by equations (13) and (l4). With these substitutions 

the value of L+ is found to be 
^n 


L t n = V (0) n + f l tf(0) n-l + f 2 w<0) n-2 * f 3 w(0) n-3 + f 0 ,,?,<0) n * 

^1 '9(0)n-l + ^2 1( ?(0 )n-2 + f*3 l(: P(0)n“3 + + 

? 2 5 n-2 + ? 3 5 n-3 + F t n + L g t (52) 

where 

F t n = e' 7te F t n . x + Jw(0) n -i + J’cpCO)^! + (53) 


L g , is one-half the tail gust force at t = ne and the f's and j's 

are defined by equation (A7) and (All), respectively, in appendix A. 


With equation (52) and difference equations (5), equations (33) 
and (34) for Lf and Mf and equation (39) for fuselage bending may 
be reduced readily to the following form: 


L fn = 7 O w (°)n + 7l w (°)n-l + 7 2 w (0) n -2 + ^3 w (°) n -3 + 7o'<P(0) n + 
7i'<P(0) n _ 1 + 7 2 ,c P(0) n _2 + 73'9(0) n _3 + 7o 5 n + h & n-l + 


7 0 5 


+ 7,8 + F + + L 


2 n-2 ' 3 ri-3 


j n 


Hr 


(54) 
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Mf n =H 0 w (°)n + ^l w (°) n -l + ^2 w (°)n-2 + ^3 w (°) n -3 + Ho ,c P(°) n + 

^1 ’9(0)11-1 + l J 2 , 9(0) n . 2 + n 3 'cp(0) n .3 + U 0 & n + % 8 n-l + 

^2 & n-2 + ^3 & n-3 ' x t F t ri - • (55) 

’•’n 

0 = r O w(0) n * r l'' ( °) n -l + r a w(0)„_ 2 ■ + r 3 v(0) n _ 3 + r 0 'q>(0) n ♦ 

ri'<p(°) n _i . r 2 'cp(0) n .2 + r3’<P(0) n _3 + r 0 S n + F 1 5 n _ 1 + 

r 2 8 n-2 + r3& n _3 + F tn + L g tn . (56) 

where the y's, p's, and r's are given by equations (A8) to (AlO) in 
appendix A. 

The complete set of loading equations as well as the fusel age 
bending equation are now available. in difference form. Equations (U9) 
and (50) apply at each spanwise station and in addition the value 
of Lf and Mf must be introduced at station 0. The coefficients t), 
v, y, and so forth are seen to involve only the physical properties of 
the airplane structure , the forward-speed and aspect-ratio parameters 
given by equations (2 6) and (5l)> certain constants derived from the 
unsteady lift function, and the time interval. The time interval e 
that is chosen should be fairly small in comparison with the natural 
period of the fundamental mode in bending of the wing. To serve as a 
guide an interval chosen near 1/30 of the estimated period of vibration 
of the fundamental mode appears to be quite satisfactory. Of course, 
some caution should be observed in the choice of this interval if the 
airplane is near a critical condition where some mode other t han the 
fundamental may predominate. For example, if the airplane is flying 
near the flutter speed, the characteristic frequency of the response may 
be near the natural torsional frequency. of the wing. The time interval 
should be modified accordingly. 

Recurrence matrix equation for response .- Equations (^9), (50), 
C5*0, and. C 55 ) for loading, equation (56) for fuselage bending, and the 
equilibrium equation (^3) niay now be combined to give the recurrence 
matrix equation for response. In order to simplify the process of 
combining these equations, only the abbreviated or symbolic form of the 
matrices which occur are used. The definitions of these matrices are 
given, unless otherwise stated, in a complete group in appendix A. 

Application of equations (^9) and (50) to each of the spanwise 
stations and of equations (5^) and (55) to station 0 leads to a set 
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of loading equations which may be put in the matrix form given by the 
following equations: 
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Equations (57) and ( 58 ) and equation ( 56 ) may now be combined to 
form the following matrix equation: 


0 


r o [ r oJ [ r o'J 


6 


r i 

r i 


rl J 


6 

1 p 1 

= 

N W M 


1 W 

+ 

l 7 i 



V 


H 

1 H 

n 

|^| [vo] M 


9 | 

n 

IH 

Vl 


VI- 


M 


? 2 

r 2_ 


_ r 2 

. 6 


r 3 - 

z 3 - 


_ r 3 *_ 

M 

12 


r »2 '] 

w 

+ 

N 

_^ 3 . 


_V 

i^l 




cp| 

n-2 

M 

V3] 


_ v 3 '_ 
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For simplicity, this equation may be abbreviated to the form: 

Ul Ul Ul Ul 



and the matrix | Lg n is defined in the section entitled "Derivation 


of Gust Forces." 

Substitution of equation (62) in equation (43) gives 

Ul b| is Ul Ul 



With the use of the notation 

00 - U - to] 

- 


( 65 ) 
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equation (64) may be written simply 

Isl 



<P 

n 


( 66 ) 


( 67 ) 


Multiplying through by the reciprocal of [d] gives finally the 
equation 



This equation gives the displacements that apply at time n in terms 
of the displacements that occurred at several preceding values of time 
(see equations ( 63 ) and ( 66 ) for the definitions of |F| n and |Q| n . 

From equation ( 68 ) the complete response of the airplane can be 
computed once the character of the gust is known. The matrix of gust- 

force values |Lgj n can be determined by the procedure given in the" ' 

section entitled "Derivation of Gust Forces." The initial conditions 
(treated in the following section) are used to obtain three successive 
initial sets of the displacements. With these sets of displacements 
the next set may be obtained by application of equation ( 68 ). With the 
newly found set and the preceding sets of displacements, the next set 
may then be found, and so forth, until a sufficient time history of 
the displacements is found from which maximum loading conditions may be 
determined. 

The reason for using the difference form of the derivatives as 
given by equations (4) and (5) might now be given. Equation (64) may 
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be considered a differential equation, since the matrix contains 

the spanwise derivative matrices CaD and. [bJ and may be likened to 
the differential equation which relates the load to the deflection for 
a beam. The unknowns are the deflections at time n. The right-hand 
terms correspond to the loading, the first term being the only one 
which is not known since it contains the unknown deflection. The sub- 
sequent inversion of the matrix [jDj leads to, in effect, the solution 
to this differential equation and, in the beam analogy, corresponds to 
the integration of the loading to obtain the deflection. When numerical 
methods are used, the deflection may be computed with good accuracy by 
integration of the loading. On the other hand, if the difference equa- 
tions which apply at an interior ordinate had been used, the matrix, [cQ 
would have appeared as an operator on one of the known deflections on 
the right-hand side of the equation. Effectively, its operation would 
be to differentiate a known deflection, corresponding in the beam 
analogy to the attempt to obtain the load which caused a given deflec- 
tion. This process, however, cannot' be done with accuracy when numeri- 
cal methods are used because of the difficulty encountered in the form 
of small differences of large numbers. The difference equations which 
apply at an outer ordinate and, consequently, lead to an integration 
process, therefore, have to be used. 

Derivation of the Initial Response 

As has been mentioned, some initial values of deflection are 
needed before equation (68) can be used. This section shows how these 
values are obtained. The airplane, just before gust penetration, is 
considered to be in level flight, and all displacements which follow 
are given relative to this level-flight condition. The initial condi- 
tions are that the vertical displacements, vertical velocity, wing 
rotation, and angular velocity are all zero. The gust force can be 
shown to start from zero and, therefore, by Newton's second law the 
additional initial condition can be established that the acceleration 
must be zero at the start of the response. These conditions can be 
satisfied, and the beginning of the response can be found by means of 
the analysis which follows. 

The initial response is assumed to be given in terms of four suc- 
cessive ordinates, denoted by w_ 2 , w_^_, w Q , and w-j_j the w Q ordinate 

is considered, as in the case of the damped oscillator, to locate the 
origin of time. The first and second derivatives at the wq ordinate 
are given by equations ( 8 ) and ( 9 ). By virtue of the initial conditions 
(the vanishing of the deflection, velocity, and accelerations at t = 0), 
the ordinate wq and the derivatives given by equations (8) and ( 9 ) 
must be zero} therefore, the ordinates w_p and w_q are found to 1 be 
related to the ordinate wq by the following relations: 
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w_2 = -8V-L 


( 69 ) 


W -1 = -w-l 


(70) 


These relations are general and must apply for deflection and rotation 
at each of the spanwise stations and for the fuselage deflection as well; 
that is, the displacements at t = -2e must be minus eight times the 
displacements at t = e, and the displacements at t = -e must he the 
negative of those at t = e . Substituting these 'conditions in equa- 
tion (64), taking n as equal to 1, and using the condition that the 
displacements are zero at t = 0 give the following matrix equation in 
terms of the displacement at t = e alone: 



(71) 


The t.em | F J-^ is zero and therefore does not appear in this equation. 

Solution of this equation gives the values of the displacements that 
occur at t = € (n = 1). 

The three sets of initial displacements required to proceed with 
equation (68) are thus known: the set of deflections found at t = e, 
the zero set at t = 0, and the set at t = -e given by equation (70), 
or simply the negative of the displacements which were found at t = e. 
In the actual case no displacements are present at t = -e, but these 
displacements may be regarded. as being of a fictitious nature the only 
purpose of which is to allow the step-by-step evaluation of the dis- 
placements to be started easily. 


Derivation of Gust Forces 

The matrix J L g which appears in the response equation (68) 

is derived as follows. From equation (15) and the notation of equa- 
tion (26), the total gust force acting over a station section at the nth 
time inverval may be given by the equation 


L g = PcZ 


0 


ne 


^ 'Kne - T )dr 


(72) 
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The integral in this expression is also of the Duhamel type and since 
the i|r function is expressed by exponential terms (see equations (21)), 
the integral may be evaluated quickly by a method similar to that 
developed in appendix D. The procedure of computing the gust force 
by this equation and then the response is not recommended, however, 
since a complete response evaluation would have to be made for each 
gust considered. Instead the procedure recommended .is to compute the 
response due to a sharp-edge gust; thenwith this response the response 
to any gust may be found directly by superposition. 

Thus for the case of a sharp -edge gust, equation (72) reduces 
simply to 

Lg n = PcZv^ n (73) 

where is the value of the i function at t = ne. This equation 

when applied to each of the spanwise stations leads directly to the 
matrix equation for gust force: 

c 0 2 0 v 0 

C 1 Z 1 V 1 
C 2^2 V 2 
c 3 Z 3 v 3 

C 4 Z 4X4 
C 5 Z 5 V 5 ; 

If the gust is uniform in the spanwise direction, the v's in this 
equation will all be equal. 

In a similar manner, one -half the gust force acting on the tail 
due to a sharp-edge gust may be shown to be 

L g tn = <Wt n (75) 

where the gust gradient is assumed to be the same as for station 0 and 

^n 

is the value of the t function for the tail. This equation and equa- 
tion (74) may now be combined to give the desired matrix |L_| as 

I sin 




follows : 
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Pt v O 0 

,0 3c 0 Z 0 v 0 

0 Pcqiq'Vq 
0 (3c 2Z2 v 2 
0 

0 pc^Zj^ 
0 Pc^Z^v^ 


*t 


(76) 


In the application of this equation it should he kept in mind 
that Lg^ does not begin to act until the tail starts to penetrate 

the gust. The time interval at which penetration starts may be taken 

x t 

as the interval nearest to the quantity — . 


Computation of Loads and Stresses 


Once the time history of the displacements has been found from 
equation ( 68 ), the normal or torque loading on the wing can be found 
with little additional work. If the notation of equation ( 66 ) is used, 
equation ( 62 ) may be written 

Ul 


n 


so 


w 


t Q 


n 


9 


n 


( 77 ) 


This equation shows that the loading matrix | P | may be found by 
adding an easily computed matrix to the matrix | Q |, the value of which 
will have already been determined in the response calculation. The 
loading matrix | P | is remembered to be defined in terms of the normal 
and torque loadings, and either of these loadings may be found 
independently of the other. 


The loadings thus found are considered to be applied statically, 
and the stresses are then found by ordinary static means. Since the 
loadings can be computed for any value of time, the complete stress 
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history of any point in the structure may he computed. In general, the 
maximum stress at different points in the structure is expected to 
occur at different times. Some- guide as to the probable time of 
occurrence of the most severe stress may be had, .however, if the 
computed wing deflection is observed. It is likely that maximum stress 
occurs in the range where wing bending appears to be most pronounced. 

The acceleration of any point in the structure may be found; if 
desired, with the aid of equation (5). 


COMPUTATIONAL PROCEDURE 


The principal results of the analysis presented in the previous 
sections are summarized herein in a step-by-step form. Only those 
steps which actually have to be performed when a determination of 
structural response for any airplane is being made are listed. In 
order to, conform with standard aircraft practice the use of inch-pound- 
second units throughout is recommended. 

The steps are as follows: 

Preliminary steps : 

(1) The wing semispan is divided into six sections and a station 

is located at the middle of each section (see fig. 3). The sections are 
proportioned in any convenient manner so that certain stations will 
fall at concentrated mass locations, such as engines or fuel tanks. 
Station 0 is located near where the fuselage intersects the wing and 
station 5 is located near the tip. The properties El, GJ, m, me, 

and mk^ are then computed at each station. 

( 2 ) From the El, GJ, and X± values determine the [aI an<j [~B ~] 

matrices by the method given in appendix C. ^ J J 

( 3 ) Compute the gust -force values at the successive time intervals 
for both the wing and the tail. (See section entitled "Derivation of 
Gust Forces.") The functions used are taken from equations (21) 

or (22) for the aspect ratios which are nearest to those of the wing and 
tail, respectively. A time interval that appears satisfactory is one 
in the neighborhood of 1/30 of the estimated natural period of the funda- 
mental bending mode of the wing. 
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The recurrence equation: 

(4) With the quantities determined in steps (l) and (2), determine 
the matrix elements given by equations (A3) to (A5) at each of the 
spanwise stations. 

(5) Compute the fuselage and tail coefficients given by equa- 
tions (A 8 ) to (All). (See definition of M-^, Mq, M 3 , M^, M<^, and M 6 

given by equations ( 35 ).) 


( 6 ) With the use of the coefficients determined in steps (4) 
and (5), set up the following matrices: £dJ , jjS-J , , [s^j , J"rJ, 

|jeJ , and |~wj. These matrices are defined in appendix A and for. the 
most part are found from simple diagonal matrices of the coefficients 
determined in steps (^) and (5). The form, for example, of the |_S^j 
matrices is illustrated in table 1 with randomly chosen numbers. All 
elements which are not shown are zero. It may be of interest to explain 
briefly the significance of the various terms that appear in the matrix. 
In order to facilitate the explanation the matrix has been partitioned 
into several submatrices. The terms in the upper left-hand box are 
associated with wing bending and the airplane vertical motion; whereas 
the terms in the lower right-hand box are associated with wing torsion 
and airplane pitching. The terms along the two subdiagonals serve to 
couple together the bending and twisting action. The terms in the 
first row and first column are associated with fuselage bending. The ' 
omission of certain terms in the matrix will lead to the matrix which 
applies to the more simple type of aircraft motion. For the case, for 
example, in which only wing bending and vertical motion are to be 
considered, computation of only the terms which make up the upper left- 
hand box is sufficient. 


(7) Determine the reciprocal of the £dJ matrix and set up the 
following matrix equation: 

Ul 
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in which 


5 


F 

-w 

F 

+ 

w 


n L J 


n-l 



In these equations the matrices containing 8, w, and cp are displace- 
ment matrices and are defined in appendix A. The matrix |F| takes 
into account the forces which develop due to the "wake effect," 
and | Lg| is the gust-force matrix which is derived in step ( 3 ). 

Equation ( 78 ) is seen to give the displacements that occur at time n 
in terms of the displacements which occurred at the times n - 1 , 
n - 2, and n - 3 . 

The initial response: • . 


(8) By use of the matrices given in step (6) and the gust forces 
which apply at n = 1, set up the following matrix equation: 



The term 



does not appear in this equation because it is zero. 


( 9 ) Solve equation (79) for the displacements. Any convenient 
method, such as the Crout method (see reference 6), may be used. The 
displacements found will be the value of displacements that apply at 
t = c or n = 1. 

The response: ' 

(10) The response may now be found by successive application of 
equation ( 78 ). The response at n = 1 has been found in step (9)j 
the response at n = 2 is next to be determined. The values of the 
displacements in the n - 2 term of the response equation are all 
taken to be zero (initial condition), and the values in the n - 3 term 
are taken as the negative of those found in step ( 9 ). The gust forces 
to use are those which apply at n = 2. The deflections that apply 

at n = 2 are then found 'by matrix algebra. For convenience the column 
matrix | Q | is evaluated first, and then multiplication of this column 
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matrix by the reciprocal of the [d] matrix gives the deflections 
at n = 2. With the newly found deflections at n = 2 and the deflec- 
tions at n = 1 and n = 0, the deflections at n = 3 can be found, 
and so forth. This process is continued until the wing bending appears 
to be the most pronounced. 

Wing loading: 

(11) With the deflections known, the value of wing loading in 
bending or in torsion can be computed directly from equation (77)- The 
stresses at any point can then be computed from the wing loading by 
ordinary static means. Since the loading may be computed at any value 
of time, the complete stress history of any point on the structure may 
be computed. 


EXAMPLE 


As an illustration of the method of analysis given in the present 
paper, the response of a typical two -engine airplane due to a sharp - 
edge gust is determined. For brevity the fuselage is assumed rigid 
and only vertical displacement and wing bending are considered. The 
weight variation over the wing semispan and the equivalent -weight 
concentrations are shown in figure 8. In this figure are shown also 
the station locations and the interval covered by each station section. 
The solution is made for a forward velocity of flight of about 210 miles 
per hour and a gust velocity of 10 feet per second. In tables 2, 3> 
and 4 are listed, respectively, the various physical characteristics 
and the factors which come from the unsteady lift function, the values 
of the i|r function and the gust-force matrix, and the matrix elements 
that are required for the solution (steps (l) to (5))> The $ function 
for an aspect ratio of 6 was chosen (see equation (20b))j and the \|f func- 
tion for an aspect ratio of infinity (equation (22)) was used. 

The [a] matrix as computed from the values of X and El 
listed in table 4 is shown in table 5(a)- In 'the computation of 
the q values shown in table 4 for station 0, the fuselage was treated 
as a concentrated wing mass. This treatment is allowable since the 
fuselage is assumed rigid and further saves the work of computing 
the 7 values (see equations (A8)). The [[A] - [Sq) ] or [d] matrix, 
which in this case applies only to bending and vertical displacement, 
is shown in table 5(h). The equation which is formed from equation ( 78 ) 
(step (7)) and which involves the reciprocal of [d] and the 

and [rJ matrices is shown in table 6. The equation for computing the 
initial response (step (8)) is shown in table 7< 



38 


NACA TN 2060 


The solution to these equations is shown in figure 9 in which 
deflection in inches is plotted against spanwise station points for 
various intervals of time. For clarity the deflections for the odd 
intervals have been left off. From these curves the consequent wing 
bending and the manner in which the airplane is swept upward by the 
gust can' be seen. The time histories of the loads (equation (77)) 
that occur at each of the spanwise stations are shown in figure 10. 
These curves indicate the presence of some second-mode excitation in 
the response. The stresses that occur at stations 0, 1, and 2 are 
shown in figure 11. The presence of second-mode excitation is, not 
readily discernible from the stress curves. 


DISCUSSION 


A method for computing the stresses and structural action of an 
airplane flying through a gust has been given. The method is based on 
aerodynamic strip theory, but over-all corrections for compressibility 
and three-dimensional effects can be made as is indicated by a suggested 
correction procedure. Some tail forces are included in the analysis 
and others might equally well be included if their character is known. 

The analysis as given is general enough to include the wing bending 
and twisting flexibilities and the fuselage flexibility. In a good 
many cases that may be considered, however, the last two of the flexi- 
bilities may prove to be of negligible importance. Some investigators 
have indicated (see reference 1) that unless the forward speed of the 
airplane approaches the flutter or divergence speed of the wing, the 
torsional deformations do not' have to be included. Likewise, in cases 
in which the fuselage is rather stiff, the effect of fuselage flexibility 
on the response may be rather small-. In such cases in which either or 
both of these flexibilities may be ignored, the analysis is^ of course, 
simplified and shortened. The example presented in the previous section 
illustrates this point. In the present state of understanding of gust- 
response analysis, enough information is not available to indicate 
definitely when and when not to include the various flexibilities of 
the aircraft structure. The analysis in the present paper may provide 
a useful means to assess their importance. The extent, for example, 
to which coupling exists between wing bending and wing torsion in any 
particular case may be seen by comparing the displacements obtained from 
the coupling terms with the displacements obtained from the noncoupling 
terms. 

Both the symmetrical and antisymmetrical types of gusts can be 
handled by the analysis given in the present paper. In general, the 
symmetrical gust is expected to produce the most severe stress condition, 
and therefore only the matrices which apply for a symmetrical case have 
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been given. These matrices were derived by using the boundary condi- 
tions for the symmetrical deformation of a free-free beam. The case 
of an antisymmetrical gust can be treated by replacing these matrices 
by the ones which apply for the antisymmetrical deformation of a free- 
free beam. The case of a general unsymmetrical gust can be handled by 
first breaking the gust into two parts, a symmetrical part and an anti- 
symmetrical part, and then treating each part independently. 

It might be of interest at this point to compare briefly the 
matrix method to a modal -function solution. One of the chief disad- 
vantages of the modal -function solution is that the modes and frequencies 
of natural vibration of the structure have to be computed in advance. 

Then, a large number of integrals which involve these modes have to be 
determined in order to set up the problem. In the present analysis 
the physical properties of the airplane are used directly in the setting- 
up of the problem. Further, in order to make the modal solution 
practical the higher modes must be dropped and only the basic or funda- 
mental modes can be used. Hence, the success of the analysis depends 
to a large degree on how. well single modal functions, one mode each for 
bending and torsion, can represent the airplane distortion. In the 
analysis of the present paper the distortions are found for all practical 
purposes as the correct values at a number of spanwise stations, at 
least to within the accuracy to which the aerodynamic and structural 
parameters are known. Also' in this analysis, probably the most 
difficult operation is the inversion of the matrix [DJ, which is 
actually not a very involved operation, especially when done by the 
quick and systematic procedure afforded by the Crout method (reference 6). 

The present paper indicates the methods for determining the 
response for both a sharp-edge gust and a gust of arbitrary shape. 

Probably the best approach, however, is to compute only the response 
■for a sharp-edge gust, since the response for any arbitrary gust may 
. thereafter be computed by means of Duhamel’s integral. To follow such 
a procedure would also save a great amount of work in the evaluation of 
the gust forces. 

One of the important features of the method of analysis presented 
is that it is not restricted to the gust problem. The approach used 
may be used to treat other problems of a similar nature. The landing 
problem can be handled by simply replacing. the distributed gust force 
by the concentrated landing forces. In the landing problem also, the 
problem is set up much more easily since the aerodynamic terms do not 
ordinarily have to be included. However, the landing problem in which 
aerodynamic forces are included may be solved by this method if desired. 
The approach used herein may also be used to solve the problem of the 
release of heavy objects such as bombs. This problem could be con- 
sidered the inverse of the gust problem; a. load is released rather than 
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encountered. Some maneuvering problems, such as the sudden deflection 
of the ailerons, and a number of other transient problems might also 
be treated by an approach similar to that giveh in the present paper. 


CONCLUDING REMARKS 


A method for computing the stresses and structural response of an 
aircraft flying through a gust has been presented. The method is based 
on aerodynamic strip theory, but compressibility and three-dimensional 
effects can be taken into account approximately by means of over-all 
corrections. The method takes into account wing bending and twisting 
deformations, fuselage deflection, vertical and pitching motion of the 
airplane, and some tail forces. A sharp-edge gust or a gust of 
arbitrary shape in the spanwise or flight directions may be treated. 

A suggested computational procedure is given to aid in the application 
of the method to any specific case. 

The possibilities of applying the method to a variety of transient 
aircraft problems, such as landing, are brought out. 


Langley Aeronautical Laboratory 

National Advisory Committee for Aeronautics 

Langley Air Force Base, Va., January 19 , 1950 
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APPENDIX A 

DEFINITIONS OF MATRICES USED IN ANALYSIS 

For convenience in presentation, most of the matrices and matrix • 
elements that are used in the analysis are defined in this appendix. 

The matrices are presented without derivations, hut their origin should 
become evident by a study of the analysis. 

Matrices.- The various matrices that are used in the analysis are 
defined as follows for the case in which the wing semispan is divided 
into six sections: (The elements which are used in the matrices are 

defined in the subsequent section.) 

w(0) 

w(l) 

v(2) 
v(3) 
w(4) 
v(5) 

9(1) 

9(2) 

9(3) 

9(4) 

9(5) 

9(6) 

8 8 

w = |w| 

9 |9| 
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k2 


p(o) 

P(l) 
P(2) 
p( 3) 
p(*0 
p(5) 

q(o) 
q(l) 
q(2) 
q ( 3) 
q(4) 
q(5) 



[*]) 

[ B ]J 


See appendix C for definitions. 
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M- 


tv] - 


M- 


M- 


^1(0) + 71 

Tli(l) 

T\±(2) 

Tli(3) 

/ 

ni(4) 
ni(5) 

hi'(O) + 71 ' 

V(1) 
n±'&) 
n ± '(3) 

' (5) 

Vi(0) + Hi 

Vi(l) 

Vi(2) 

V ± (3) 

Vi (4) 

v i(5jj 

Vi'(O) + Hi' 

Vi'(l) 

Vi' (2) 

Vi' (3) 
v ± ' (4) 

v i'(5)| 
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LJJ = [j 0 0 0 0 0 j 


[j'J =[_j'00000| 
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j [jj [yj 

0 [ g ] [g] 


Matrix elements.- The matrix elements which appear in the matrices 
defined in the previous section are expressed foy convenience in terms 
of the following common factors: 

\ 


3 = m A *pU 



®0 = a l 
*0 = -^ a l 


(Al) 


in which the last four are associated with the 4> function for the 
wing. (See equation (23).) With these factors the elements that must 
be computed at each spanwise station are as follows: 


do = -<&o)PcZ + pc I (®Q + 

d x = |(1 - $ 0 )pc Z 


*2 '= -^(1 - *o)PcZ 
d^ = ^^1 - $o)PcZ 

dQ 1 = 57(1 - ®0)Pc 2 l(| - |) + pcZ £(1 - * 0 ) 

§®0 U€ - {*0 + h*) C (l - f)] 



dl * =-|(l - $o)Pc 2 z(3 _ 

d 2 - = ^(1 - %)Pc 2 z(| - 
d 3' = "Ti (l ' ^ )pc2i (t _ |) 


(A2) 
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The coefficients which must he computed for the fuselage and tail 
are expressed in part in terms of the following common factors: 


Pt = 2 ffl A t ItpUS t 


y _ 2U . 

" cSt Xt 


$ to = 


% - -?t a l t 


®to - ?t^ a lt J 

in which the last four are associated with the 4 functions appropriate 
to the tail. Also used are the generalized masses given by equa- 
tions (35) and the value 0]_ as given in equation (l4). With these 
factors’ the coefficients for the tail and fuselage are as follows: 


f 0 =-£(1 - »to)et + Pt(®to + |V) 

f i - K 1 - »to)H 
f 2 --U 1 - **oK 
f 3 - U 1 - \Y t 


- %)(*t ♦ ic t )s t ♦ 0t[u(l - %) - 
5%“' - (*to + 5*tO S )( x t + rt)] + IliPtct 
-K 1 • 1 ' t o)( xt + - n® tc t 

A( X ‘ + + 

" 3^( 1 ‘ *tj ( x t + £t)st - 


(A7a) 


(ATb) 



50 


MCA TN 2060 


% - %) (1 + gct?i)pt - Ptfei(i - *to) - 

i$ 0 u "l " (®to + |^to 6 )(*- + |°t e l)] " ife? tC t e l 
fl = l(l _ 3^) (l + ict 0 l)pt + T^Ptct^l 
f 2 = -^(l - ^toX 1 + ^t0l)pt - ^Ptctei 
^3 = J^ 1 " + |°t e l)pt + jJ^tctPl 


2Mi 

7o = -— + f 0 

e d 


5Mi 

£ 

lain 


r l ■ -t + f l 


7 2 - - — + f 2 

Mi 

r 3 - i + f 3 


V 

7i' 

V 

7 l' 


2Mo 

i# +f o’ 


5M 2 

5 +f I* 


6 ' 

7T 


+ tr 


Mo 

■ ~ + f 3 ' 


> (A8a) 

y 


> (A8b) 


J 
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2M 3 _ 

-5M,- 

r i - w + fl 

ItM 

^=- e 2 
Mo 

7 3 = + f 3 


3,; 
+ f 2 


Ml = - x t fl 


2Mo 

mo = -%r ~ x t f o 
5m 2 

r 

6 

l4Mo 

M2 = — - x t f 2 
6^ 

M2 

^3 " -3 “ x t* 3 


2 M 4 

M0 ' = ~^2~ ' x tV 

5Mj^ 

mi' = -gr - x t fi' 
4fJl 

M2' - x t f 2* 

M 4 

^3’ = ^2 - x tf 3 ’ 
2M 5 

Mo = -72“ " x t f o 

Ml^-^-X^ 

€ d 

M2 = ■ x t f 2 

M 5 _ . 

M3 = --2 - x tf 3 


> (A 8 c) 

y 

\ 

> (A9a) 

J 

> (A9b) 

/ 

> (A9c) 

y 


/ 
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2Mo 

r ° = + f ° 


5M^ 

r i = 7? + f i 




r 0 = 


+ fr 


r, = ^2. + f_ 

3 € 2 3 


2Mc 

»' - -£■ + f o’ 


. = _^ + 


fi* 




r 2’ " -f + f 2' 

Me: 

r 3 '--^ +f 3 ' 


2 Mr 5 

r o = ~~gr + - “f m 6 

5 Mg _ 

r-L = -=2. + f 
? 2 = 


e 2 - 1 
— 3 “ + ? 2 


M 6 _ 

r 3 = * + f 3 

- 7 t€ 


j - Mt 0 ee 

y = pt ee ”^ te jj- 


u$ to - $ to(x t + -C t 


*\ 


\ 


y 

\ 


i 


$ 


i = Pt 66 


- 7 1 € 


[>Vi + %(i + ^i)] 


y 
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(AlOa) 


(AlOb) 


(AlOc) 


(All) 
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APPENDIX B 

DERIVATION OF DIFFERENCE EQUATIONS 


In this appendix the parabolic and cubic difference equations for 
the first and second derivatives of a function are derived. 

Parabolic equations .- For the parabolic difference equation, 
consider the function shown in figure 12(a). This function is assumed to 
be replaced by the arc of a parabola which passes through the three 
ordinates a, b, and c. It can be verified readily that such a curve 
can be given by the equation 

“ - f(f - 1 )(? - 2) -’$-*)♦!#- 1) w 

The first and second derivatives of this equation at y = e are given 
by the equations 


dw l m c -a 

dyjy=e 2e 


(B2) 



2b + a 
e 2 


(B3) 


These equations are the standard difference equations for the first 
and second derivatives of a function. The derivatives are purposely 
taken at the middle of the three ordinates because the resulting 
equations are suitable for use in the simplification of many problems. 

If the derivative had been taken at an outer ordinate, the approximation 
afforded would not be accurate enough to be useful . 

Cubic equations .- The cubic difference equations may be derived 
in a manner similar to that for the parabolic equations . In this case 
four successive ordinates are used. (See fig. 12(b).) The function 
is replaced by a third-degree curve which is given by the equation 
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5b 


“ - - i(? - i)(f • *)(? - 



+ 



W 


Because of the increase in accuracy that results from the use of a 
higher-degree curve, the first and second derivatives may he taken at 
an outer ordinate with an accuracy which is about equivalent to that 
given by equations (B2) and (B3) . The derivatives at y = 3e are 


dw lid - 18c + 9b - 2a 

dy y=3e " - 


(B5) 


d^w 

dy 2 

_ 


y=3e 


_ 2d - 5c + j£b - a 

6 2 . 


(B6) 


If taken at the third of the four ordinates, the derivatives are 


dw 


y=2e 


_ 2d + 3c - 6b + a 

_ — 5 - 


(B7) 


d^w _ d - 2c + b 

iy^|y=2e ? 


(b8) 


Equations (B 5 ) and (b 6) are used in the derivation of the response, 
equation for an airplane in a gust. Equations (B 7 ) and (b 8) are useful 
in the derivation of the initial response. 
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APPENDIX C 

DERIVATION OF MATRIX EQUATIONS OF EQUILIBRIUM 


In this appendix the matrix equations 


H 

|w 

-m 

[ B ] 

hi 

1 ■ M 


(Cl) 

(C2) 


for symmetrical bending and twisting of a free -free beam under normal 
and torsional loads are derived. 


Bending . - In accordance with the assumptions made in this paper the 
wing semispan is considered to be divided into six sections with a station 
point at the center of each section (see fig. 3). The inertia force of 
the mass and the aerodynamic force that develops over each section is in. 
turn assumed to be concentrated at the respective station points. The 
wing is thus effectively a beam bending under six concentrated loads and, 
as such, will have a linearly varying moment between each station. The 
following general equation for the moment between the i and i + 1 station 
may therefore be written: 


M = a^ + b^y 


where 


a i = 


H' 


1+ ^7rl M < 1) + 11 


b i-bxI7l[ M(1 +1) 


(C3) 


in which y(i) is the abscissa to the ' i station. 


The wing is further assumed to have a linear l/EI variation 
between stations with the correct value of l/EI at each station. This 
type of variation would lead to an El curve which follows very closely 
the true stiffness curve of the wing and which of course has the correct 
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values of El at each station. A general equation for l/EI may 
therefore also he written; thus, 


5 ' 01 + iiy 


(C4) 


where 


c i 



1 

El(i) 


y£i) 1 

bX i+l El(i + 1) 


d 1 = 


1 

1 

"k^i+1 

El(i + 1) 



With equation (C3) and equation (C4) the well-known expression 
relating deflection to moment for a heam may he written 


= (*i + hj.y)(ci + dj.y) 

dy 2 EI 


(C5) 


The deflection may he found most conveniently . from this equation hy use 
of the engineering team theorem which states that the deflection of one 
point on a beam relative to the tangent of the deflection curve at 
another point is equal to the moment about the displaced point of the 
M/EI diagram between the two points. ' In this case symmetrical loading 
is being considered and therefore the boundary condition at the center 
line is that the slope must be zero; the deflection of each station 
relative to this point therefore may be readily computed. Fortunately, 
because of the convenient analytical representation of M/EI, these 
deflections may be found by exact integration. The deflection, for 
example, at station 4 due to the M/EI variation between stations i 
and i + 1 may be given by the expression: 


py( i+1) 
'y(i) 


(aq + biy)(ci + d ± y ) |y ( 4 ) - y]dy 
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Consideration of all the expressions of this sort leads to the total 
deflection of each station relative to the wing center line. From this 
deflection the more useful deflection relative to station 0 can be 
readily determined. The values of the deflection thus obtained are 
found to be expressible by the following matrix equation: 


v(l) 


w(2) 


v(3) 

w(4) 

_ b2 

EI(0) 

w(5) 



a ll 

a 12 

0 

0 

0 

®21 

®22 

823 

0 

0 

a 31 

a 32 

a 33 

a 34 

0 

a 4l 

a 42 

a 43 

a 44 

a h5 

a 5l 

a 52 

a 53 

a 54 

a 55 


M(0) 

M(l) 

M(2) 

M(3) 

M(4) 


(C6) 


where the matrix elements are defined by the equations: 


p 

^*11 = ^0^1 ^*1^1 

&P2. = * 0 (*1 + Xp) + Xj^A^ * 

a 31 = + ^2 + ^3) + ^i^Ai + + >^)% 

= * 0 (*1 + ^2 + ^3 + + + + X^ 

a 51 = + ^2 + ^3 + + ^5) + + ^(Xg 


ACTa) 


+ X^)B^ 

+ x^ + x^ + ^5)^1 


a l2 = X 1 2 C 1 
a 22 = XqSCi 
a 32 = X 1 2 C 1 
2 

a^p = X]_ C-l 

a 52 = X^C^ 


p 

+ X^XpD-^ + Xp Ap 

+ ^l(^P + ^3^1 + Xp^A p + XpX^Bp 

+ X]_(Xp + X^ + + Xp^Ap + Xp(X^ + X||)Bp 

+ X]_(Xp + X3 + Xl| + X^)D]_ + Xp^Ap + Xp(X^ + Xlj. + Xcj)Bp 


>(C7b) 
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a 2 3 — ^2 ^2 

a 33 x 2 ^2 + X 2 X 3®2 + ^3^3 

a 43 = ^2 ^2 + x 2^ x 3 + H) D 2 + ^3^3 + ^'•^'4®3 
♦ 

a 53 = ^2 ^2 + M x 3 + x 4 + x ^)^2 + x 3 2a 3 + ^3(^4 + ^5)^3 


(CTc) 


9/^1^ — X^ C-^ 

a 44 = x 3 2 ^3 + X 3 X 4®3 + ^4 2 A^ 

a 54 = x 3 2c 3 + ^3(^4 + ^5)^3 + ^4 2 A^ + 


(C7cL) 


a 45 = H 2 °4 


a 55 ~ x 4 2 °4 + x 4 x 5 d 4 + X 5^ A 5 


(C7e) 


in which 


A, 1 

+ i 1(0) 

f 1(1 - 1) 

12 5 ( 17 . 

B- - 1 K. 0 ) 

1 1(0) 

■“I “ 3 

1(1 - 1) 

6.1(1) 

c , - i- I<0) 

. + 1 I (°) 

°i - To 

2 I(i - 1) 

12 l(i) 

r, _ 1 1(0) 

. 1 1(0) 

1 « Ki-D 

3 'I(i) 


r (c8) 



The moment M(5) does not appear in equation (C 6 ) because the boundary condition is used that the 
moment at station 5 is zero. For convenience equation (C6) may be given in the abbreviated form: 
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Substitution of equation (Cll) into equation (C9) gives 



(C12) 


Multiplication through, by the reciprocal of 
the equation: 


b 3 

El ( 0 ) 



results in 



(C13) 


This equation thus gives the loads in terms of the deflection of each 
station relative to station 0. In the case under consideration^ however 
the wing is a free body capable of motion through space and therefore to* 
set up properly the equations of motion the deflection must be given 
relative to a fixed datum line. This datum line is most conveniently 
located as the position of the wing prior to action of the disturbing 
loads. Consideration of the following sketch 



will show therefore that the following relation must exist: 

w = w - w(0) (CliO 

Substitution of this equation into equation (C13) gives 



(C15) 
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6l 


El(O) 


To aid in the derivation — [MM] 
general form: 


-1 


is now written in the 


El(0) 

v3 


[MM] 


-1 


b n 

b 12 

*13 

*14 

*15 

*12 

h 22 

*23 

V 

*25 

h i3 

*23 

*33 

*34 

*35 

*14 

*24 

b 34 

*44 

*45 

*15 

*25 

b 35 

*45 

*55 


(C16) 


Thus with this equation, equation (C15) may he transformed to the form; 


where 


’01 

*11 

*12 

*13 

*14 

*15 

’02 

*12 

*22 

*23 

*24 

*25 

’03 

*13 

*23 

*33 

*34 

*35 

*04 

*14 

*24 

*34 

*44 

*45 

*05 

*15 

*25 

*35 

*45 

*55 



w(0) 



w(l) 


p(l) 


v( 2) 


P(2) 


w(3) 

= 

P(3) 


w(4) 


p(4) 


w(5) 


p(5) 


b 01 = "C b H + b 12 + b 13 + *14 + *15) 


3 02 


3 03 


-(h 


"^12 + ^22 + ^23 + "^24 + 


"b-^3 + D 23 + 1) 


>5) 

33 + h 3 4 + h 3 J) 


> 


b o4 - 
*05 = 


+ b 2h + b 34 + b 44 + b 4^) 
■(*15 + *25 + *35 + b 45 + b 5^) 


(C17) 


(Cl8) 
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Equation (C17) is noted to. express all the loads except p(0) in 
terms of the six deflections. An additional equation in which p(0) is 
expressed also in terms of the deflections may he established by use of 
the condition that all the loads acting on the wing semispan must add 
up to equal zero; that is. 


p(o) + p(l) + p(2) + p(3) + p(4) + p(5) = 0 (C19) 


This condition automatically satisfies the two boundary conditions that 
the shear must be zero at the tip and center line of the wing. Thus if 
the five equations represented by equation (C17) are added, and use is 
made of equations (Cl8) and (C19), the following equation results: 


t 0 o w (°) + t> 01 w(l) + ^ 02 w( 2) + b 03 w(3) + h 0 4w(4) + b 05 w(5) = p(0) (C20) 


where 


J 00 




01 


+ *02 + * AO + * rt1, + * 


03 " 04 


)s) 


(C21) 


This equation may now be combined with equation (C17) to give finally 


b 00 

*01 

*02 

*03 

*04 

*05 


w(0) 


p(o) 

*01 

b n 

b 12 

*13 

*14 

*15 


w(l) 


P(D 

*02 

*12 

*22 

b 23 

*24 

b 25 


v(2) 


P(2) 

*03 

*13 

^23 

b 33 

V 

*35 


w(3) 


P(3) 

*04 

*l4 

*24 

*34 

*44 

*45 


w(4) 


p(M 

*05 

*15 

*25 

b 35 

b 45 

*55_ 


w( 5) 


p(5) 


(C22) 


This equation is thus the desired matrix equation which relates the normal 
loads to the deflection. If the square matrix is denoted by M , the 
equation may be abbreviated conveniently to the form 


[a]M=|p| 

which is the form used in the text. (See equation (40.). 


(C23) 
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As an aid in computational work, a summary of the steps involved in 
the determination of is given to close this section: 

(l) From the I values at the respective stations compute the 
coefficients given ty equations (C8) 


(2) With these coefficients determine the matrix elements given- ty 
equations (C7). These elements form the matrix which is defined 

by equations (Co) and (C9) 


(3) Multiply the [hJ matrix by the matrix, which is defined 

by equations (CIO) and (Cll). The result should be a symmetrical matrix; 
this property serves as a very useful computational check. 


( k ) Invert the 


b3 


El(0) 



% H 2 matrix. This matrix should also be 


\ ~ / l_T — — — * 

symmetrical. (The Crout method (reference 6) serves as a rather quick 
and useful means for performing the inversion. ) 


(5) Add the columns of the inverted matrix and place the negative 
of these sums at the top of their respective columns such as to form a 
new' row of matrix elements. Then add these sums and place the negative 
of the sum as the first matrix element of the newly formed row. A new 
column headed by this value is thus in the making. Fill in the remainder 
of the column with the respective elements of the new row; that is, the 
appropriate values should be inserted to make the matrix sy mme trical 
This final matrix is the desired matrix . 


Torsion . - For the torsional case the torque loads q are assumed 
to be concentrated at the stations just as in the case for the normal 
loads p. Consideration then of the following example torque -diagram 
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will show that the following equations must apply: 


q(0) = - T(l) 
q(l) = T(l) - 1 ( 2 ) 
q(2) = T(2) - T(3) 
l(3) = T(3) - T(4) 
q(4) = T(4) - T(5) 
q(5) = T(5) 


(C24) 


where T(i) represents the total torque present in the i interval. No 
torque exists between the wing center line and station 0. 

To aid in the derivation, the assumption is made that l/GJ varies 
linearly between stations. A typical T/GJ diagram between, say, the 
i - 1 and the i station would appear as follows: 


T(i) 

Gjfi - 1) 

i - 1 i 



jtl) 

Gjfi) 


From the differential relation |3> = jL the fact may be observed that 

ay GJ * 

the change in angle of twist between two stations is equal to the area 
of the T/GJ diagram between the two stations; therefore, 


bX.j 

T(i) 

2 

GJ(i - 1) 


T(i) 

GJ(i) 


(C25) 
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If the notation 


Ji " 


2G 

bXi 


(C26) 


J(i - 1). J(i) 

is employed, equation (C25) nay he written 


T(i) = jiji(i) - qp(i - D] 


(C27) 


Application of this equation to each of the spanwise stations gives the 


T: 

- 

t(i) = ^[9(1) - 9(0)] 
T(2) = J 2 [<P( 2 ) ‘ 9 (!)] 

T( 3 ) = J3 

[cp( 3 ) - 9(2)J 

T( 4 ) = 

| 9 ( 4 ) - 9 ( 3 )J 

T( 5 ) = J 5 

[ 9 ( 5 ) - q( 4 )J 


(C28) 


Substitution now of these equations into equations (C24) gives the 
desired equations relating the torque loads to the angle of twist. The 
equations thus found can be given in the matrix form: 


J 1 

~h 

0 

0 

0 

0 


9(0) 


4(0)' 

"Jl 

(Ji + ^2) 

-J2 

0 

0 

0 


9(1) 


q(D 

0 

~h 

(J2 + J3) 

"J3 

0 

0 


9(2) 


1(2) 

0 

0 

" J 3. 

O3 + Ji|) 

-J4 

0 


9(3) 


l(3) 

0 

0 

0 


( j h. + J5) 

~h 


9(4) 


4(4) 

0 

0 

0 

0 

~b 

h 


9(5) 


l(5) 


(C29) 
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which can be abbreviated to 


[b]M = |q| (C30) 

the form used in the test. (See equation (4l).) Thus all that is 
involved in the computation of the matrix (Bj is the evaluation of the 
matrix elements by means of equation (C26). 
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APPENDIX D 

RECURRENCE EQUATION FOR THE EVALUATION OF DUHAMEL'S INTEGRAL 
INVOLVING AN EXPONENTIAL KERNEL 


The derivation of ‘a rather simple recurrence relation for the step- 
by-step evaluation of the three unsteady lift integrals appearing in 
equation (25) is presented. This derivation is made possible because 
the kernels of the integrals are expressible in exponential form. 


From equation (23) the first and second derivative of the 
G> function may be written 


<t» = 


-A^t 

2U c o . 

- — Xa, e = 0 © 

c 0 1 0 


■7t 


(EL) 


$ = 


oo 2 1 


.. -rt 

V 


(E2) 


where 


r = xS 

c o 


»o • ' ra i 


*0 " 


With these equations the three integrals of equation (25) may be 
combined conveniently into the following single integral denoted by 1^ : 




cple 


-7(t- T ) 


dr (D3) 
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For convenience the notation 


Y -= |$ 0 PcIw - <|> o (3cIU + $ 0 Pc 2 Z^ - ? )]J 


is introduced and thus equation (D3) becomes 


■/ 

J 0 


t -r(t-T) 
Ye 


cLt 


0 *) 


or 



(D5) 


Mathematically, the integral in this equation may be interpreted to 
represent the area under the function given as a product of Y and e ? T . 
In accordance with numerical evaluation processes, the interval 0 
to t may be divided into a number of time stations of interval e . 

The product of Y and e 7T may then be found at each of the time 
stations and from these products the area under the curve may be deter- 
mined in first approximation by the trapezoidal method of determining 
areas. Thus, if the n time station corresponds to time t, the 
expression for 1^. may be approximated as follows: 


~ “ ee 


-yne 


ye 


y2e 


V 


+ Y 2 e 


Y e 
n-1 


y(n-l)e 




Tne 


n 


(D6) 


where Y Q does not appear since the initial conditions are. used that 

the deflection w and rotation cp are zero at t = 0, and' therefore 
Yq is zero. (See equation (DU).) More accurate methods, such as 

Simpson's method, could be used for determining the area under the curve, 
but because of the small interval chosen the consequent increase in 
accuracy is negligible. If the notation 


F n = ee 


■/ne 


Y^ 76 . Y 2 e 72£ ..... Y n . 1 e 7(n_1)6 


(D7) 


is introduced, equation (d 6) may be written Bimply 


I 


n 


= F n 



m 
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If equation (D5) is expanded similarly, only for an upper limit 
of t - e, the expanded result would he 


'/(n-l)d 


I = ee 
n-1 


„ 7e y2e 

fY^e + Y 2 e + 


+ Y 
n' 


y(n-2)e 1 r(n-l)ej 

■2 6 + 2 n-l® 


(D9) 


By analogy with equation (D7), however, 


F n-1 = ee 


-7(n-l) € 


+ Y 2 e 72e + •** + Y n -2 e7(n " 2)e (DIO) 


and therefore equation (D9) Becomes 


T n-1 " F n-1 + 1 Y n-1 (Dll) 

A study of equations (D7) and (DIO) shows that the following relation 
must exist: 


F n = e F n-1 + ee Y n-1 ‘(D12) 

Now^ if equation (D4) is used to rewrite and Y . in equations (D8) 

and (D12), the value of may he given finally hy the equation: 


where 


n 


= F. 


n 



V cZeW n 


- |pcZe 





(D13) 


F n = e ’ 76F 


-7c 


'n-1 + <5>o ee Pc2w n-l ' 


-ye 


U0„ + c 


(? - f )*0 


^n-1 


(Dll) 


The value of the unsteady lift integrals is thus given hy 
equation (D13). As regards the analysis given in the present paper, 
"n-1 8111 9n-l 8X6 the values of deflection and rotation which have, 

sa Y> just heen detennined from the recurrence equation for response. 
The value F n _-]_ was also established and therefore F n can he deter- 
mined as a definite quantity. The value I n is thus seen to he given 



7 ° 

in terms of the known F n and in terms of w n and cp n 
next values to he evaluated from the recurrence equation 
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APPENDIX E 
MATRIX ALGESIA. 

This appendix is written for "those not familiar with matrix notations 
or matrix methods. All the matrix algebra necessary for the understanding 
of this paper is described hereinafter by way of examples. 

Matrix definition . - Some of the basic types of matrices are illustrated . 
by the following arbitrary matrices which are of the third order: 

The column matrix 

2 
1 
-1 

The row matrix 

[2 -3 xj 

The square matrix 

2-3 1 

1 2-2 

-1 -1 3 

The diagonal matrix 

k 0 0 

0 ' 3 0 

0 0-1 

The identity matrix 

1 0 0 

0 1 O' 

0 O' 1 
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Element definition . - Each of the terms that appear in a matrix is 
d.©fine4 as an element. Its position is usually denoted, in a row by the 
number of terms from the left and in a column by the number of terms 
from the top. 

v Matrix addition . - The addition of two matrices produces a single 
matrix. Addition is performed by simply adding together corresponding 
elements. For example. 


2-31 


1 

0 

rH 

1 

-3- 

1 


6 -4 1 

1 ' 2 -2 

+ 

0 3 2 

= 

15 0 

_-l -1 3_ 


Vo 

0 

1 

H 


1-12 


Multipli cation of a matrix by a scalar number . - In the multiplication 
of a matrix by a scalar number every element in the matrix is multiplied 
by the number. For example, 


2-3 1 


4-6 2 

12-2 

It 

2 4-4 

: l -1 3_ 


1 

VO 

CVJ 

1 

CVJ 

1 J 


Multipli cation of a column matrix by a row matrix .- The product of 
a column matrix and a row matrix is equal to the sum of the products of 
the correspo n d in g elements . For example, 


L 2 


-3 


ij 


2 

1 

-4 


= ( 2 x 2 ) + (-3 x 1 ) + [l + (-4)] = -3 


Multiplication of a column matrix by a square matrix . - The multiplication 
of a column matrix by a square matrix produces a column matrix. Consider 
the following set of three simultaneous equations: 


2y x - 3y 2 + y 3 = 

y l + 2y 2 " 2y 3 = a 2 V 


- y-, 


y 0 + 3y 0 = a 


3 


(El) 
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The procedure adopted in matrix algebra is to write these equations in 
the matrix form 


2-31 


7 1 


a l 

cu 

1 

OJ 

H 


CM 

= 

a 2 

[-1 -1 . 3 


7 3' 


a 3 


(E2) 


where the multiplication of the |y| column matrix by each row in. the 
square matrix produces the respective elements in the |a| column matrix. 
(See multiplication of a column matrix by a row matrix.) 

In order to simplify the presentation of an analysis, the symbolic 
or abbreviated matrix form is used quite often. The symbolic form of 
equation (E2) is simply 


MM = M (E 3 ) 


The determination of |a| by the multiplication of |y| by [m] is 
illustrated with arbitrary values of y, say y^ = 4, y^ = 5, and y = 6, 
by the following equation: 


2 -3 

1 2 

-1 -1 


( 2 X 

( 1 X 

(-1 x 


4) 

4) 

4) 


+ (-3 x 5) 
+ (2x5) 
+ (-1x5) 


+ (1x6) 
+ (-2x6) 
+ (3x6) 


-l 


a 



1 

2 

= 

*2 

9 




(E4) 


Multiplication of a square matrix by a square matrix , - The multi- 
plication of two 3quare matrices produces a square matrix. Multiplication 
is performed by letting the multiplying matrix operate, as in the' preceding 
section, on each of the successive columns in the matrix being multiplied 
to produce corresponding successive columns in the product matrix. For 
example. 


'2 -3 1 " 


1 _2 3 ~ 


-2 -12 14 " 

1 2 -2 


13-2 

= 

5 2-5 

-1 -1 3 


-11 2 


-5 2 5 


(E 5 ) 
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Order of multiplication .- In general the commutative multiplication 
law of ordinary algebra does not hold in matrix methods; that is. 


a| |b I 4 |b| |a 


Therefore, whenever the product of several matrices is indicated, these 
matrices must be multiplied together without changing their order. 

1 

Matrix partitioning and submatrices .- A matrix may be partitioned 
or divided at will into smaller matrices.. For example, the left- hand 
side of equation (E4) may be partitioned as follows: 


2 | -3 1 


4 

1 [ 2 -2 


5 

l 

-li-l 3 


6 


The matrices which are formed by the dividing lines are called submatrices. 
These submatrices may be treated as though they were elements when matrix 
operations are performed. For example, with the notation 



the multiplication of the foregoing partitioned matrix is as follows: 


2 ' a 

1 . 


4 


8 

+ ad 

l 

b 1 c 

1 


d 


4b 

+ cd 
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The reciprocal of a matrix and the Identity matrix . - ordinary 
algebraic methods the formal operation involved in the solution for x of 
the equation 


mx = a 

is the multiplication through by the reciprocal of m; thus, 

x = m -d a 

The same ( f ormal operation may be applied to matrix equations. For example, 
the solution for |y| • in equation (E3) is simply 


kl = M _1 

where [m] is the reciprocal, or the inverse, of [mJ . 

- The reciprocal of a matrix is found as the matrix which satisfies 
either of the equivalent equations 

[m ]- 1 [m] - [l] 

[m] [m] -1 = [1] 

where [ij -is the identity matrix. For equations (E2) and (E3), the 
reciprocal of [M] is found as the matrix which satisfies the equation 


"2 -3 1 


b l C 1 d l 


10 0 

1 2 -2 


b 2 c 2 d 2 

II 

0 10 

-1 -1 3 


b 3 c 3 d 3_ 


0 0 1 


If this equation is considered in relation to equations (El), (E2), and (E5), 
the values b^, bg, and b^ would simply be values of y^, yg, and y^ which 

satisfy equation (El) for a-|_ = 1, a 2 = 0, and a^ = Oj c^, c 2 , and C3 would 

be the values for a-^ = 0, a 2 = 1, and a^ = 0; and dj_, dg, and d^ would 

be the values for • a^ = 0, ag = 0, and a^ = 1. For this example, the 

solutions are 


’ 1= 3 

Cl = — 

1 3 

d l 

_ ± 
3 

'2 = 

c p -i 
2 12 

*2 

Ha 

II 

, _ 1 
3 12 

n = -L 

d_ 


3 12 

3 

12 
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The Crout method (reference 6) provides a very quick and convenient 
method for determining these solutions . 


The determination of y by the operation [m] on |a| is illustrated 
as follows for a^ = -1, a 2 = 2, and a^ = 9 : 


1_ 

12 


' k 8 k 


-1 


4 

-17 5 


2 

= 

5 

1 

O- 

lT\ 

H 

1 


9 


6 


(Ell) 


The operation performed hy this equation can be seen to be the inverse 
operation of equation (E4) . 
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TABLE 1 

ILLUSTRATION OF THE [s] MATRICES 
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TABLE. 2 


PHYSICAL CHARACTERISTICS AND UNSTEADY LIFT 
FACTORS FOR EXAMPLE AIRPLANE 


b, in , 560 

c o> 154- 

EIq . ' 

lb/in. ' .165 

p, lb/ft 3 . 0.0765 

fmph 210 

u > [in. /sec 3700 

v, in. /sec 120 

e, sec 0.01 

As, half-chords ....... 0.48052 

M 0.276 

A 10 

m A 0.861 

P . . . 0.0011^7 

0.381 

7 .18.3078 

e - ? 6 ....... 0.832703 

$ 0 = a i • • 0.361 

®o -6.60912 


"0 


120.9983 



TABLE 3 


t ORDINATES AND GUST -FORCE MATRIX FOR 
EXAMPLE AIRPLANE 


n 

(Equation (22)) 

0 

0 

1 

.22105 

2 

.36744 

3 

.46716 

k 

.53741 

5 

.58888 

6 

.62830 

7 

.65980 

8 

.68594 

9 

•7o84o 

10 

.72820 

11 

.7^595 

12 

.76215- 

13 ' 

.77707 

14 

.79086 

15 

.80373 

' 16 

.8157^ 

17 

.82696 

18 

.83748 


J g 


= 120 


17.81*04 

15.7552 

12.1811 

10.5295 

8.77455 

7.01964 


k 


n 
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w 

>4 
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00 rH 

rH 

O 

P“ 
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rH 
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rH OO rH t VO uN 

t*— O O VO 0-4 

ON -4 0 ON Q CO 
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O LPs LTN 1 — 1 rH VO 

• •*••• 
VO H t— CVJ rH 1 

lf\ on 1 1 1 

1 1 

B 

H H O O ' O O 

O O On On On On 

rH 1 — 1 

0 

-4 VO CO CVJ LTN CO 

IT\ on rH O CO VO 

rH rH rH rH 

EI 

LTN 

0 0 0 0 m\ n— 

0 n- 0 cvj oj vo 

t — CO rH n OJ 

CVJ H rH 

n 

ON CO n- VO VO VO 

O rH rH 1 — 1 rH rH 

«••••• 

O 

Station 

0 H OJ on -4 LTN 
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. TABLE 5 

[d] matrices for example airplane 

(a) The [a]| Matrix. 
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on m 00 


00 rH 

H 00 CV1 VO CO 

cvj on n- cvj vo 

co 4 on 4 


CVJ On 00 t— VO 
I CD 4 00 
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in 
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on 
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co 


in 
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On 
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On 
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ON 
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-4 

rH 

on 

on 

d 

ON 
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rH 

co 

co 

on 
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-4 
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vo 

on 
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0 

CVJ 

in 

in 
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1 

•N 
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RECURRENCE EQUATION FOR RESPONSE OF EXAMPLE AIRPLANE 
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TABLE 7 
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Figure l . - Damped oscillator and suddenly 
applied force . 



O .05 JO J5 .20 

t\ sec 


Figure 2.- Comparison of exact and difference- 
equation solutions for response of damped 
oscillator. 
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Figure 3.- Division of wing into sections. 



figure 4. - Displacements of a wing cross section. 



Figure 5. ~ Coordinate system for fuselage displacement. 
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Figure 6.- Lift functions for sudden change in 
angle of attack. (See equation (20).) 



Figure 7.- Lift functions for wings entering a 
sharp - edge gust. (See equations (2!) and (22). ) 
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for example two - engine aircraft. 



Figure 9.- Response of example airplane due to / O-foot - 
pet second sharp - edge gust. U= 210 miles per hour. 
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Figure 10. -Time history of station loads for 
example airplane. 



Figure //.- Bending stress developed in example 
airplane due to 10' foot-per- second sharp - edae 
gust. 
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( o) Para bolic . 



Figure /2r Functional notation used in the derivation 
of parabolic and cubic difference equations. 
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